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Discretizations of the Feynman-Kac path integral representation of the quantum mechanical density 
matrix are investigated. Each infinite-dimensional path integral is approximated by a Riemann 
integral over a finite-dimensional function space, by restricting the integration to a subspace of all 
admissible paths. Using this process, a wide class of methods can be derived, with each method 
corresponding to a different choice for the approximating subspace. The traditional "short-time" 
approximation and "Fourier discretization" can be recovered from this approach, using linear and 
spectral basis functions respectively. As an illustration, a novel method is formulated using cubic 
elements and is shown to have improved convergence properties when applied to a simple model 
problem. 
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The path integral approach provides a powerful 
method for studying properties of quantum many-body 
systems [Q. When applied to statistical mechanics 
each element of the quantum density matrix is expressed 
as an integral over all curves connecting two configura- 
tions: 
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The symbol T> [x (t)] indicates that the integration is 
performed over the set of all differentiable curves, x : 
[0, (3h] R'^, with x (0) = a and x (/37i) = b. The integer 
d reflects the dimensionality, with d ~ 3N for a system 
of iV-particles in 3-dimensional space. The functional $ 
can be derived from the classical action by introducing 
a relationship between temperature and imaginary time 
{it — (3h) 1^]. In this Letter, we will restrict our attention 
to the quantum many-body system, for which $ takes the 
following form: 



<i>[x(r);/3]= / -J2 

•^0 ^ ^=l 



m, ±, (r)2 + [x (r)] dr. (2) 



Calculating the path integral in (|i|) is a challenging 
task, which in general cannot be performed analytically. 
It is only for simple model problems, such as quadratic 
potentials that an exact solution can be obtained. For 
more challenging systems, the path integral has tradi- 
tionally been estimated using either the "short-time" ap- 
proximation (STA) or "Fourier discretization" (FD) 
1^,^. Many authors have proposed improvements to the 
standard STA and FD, using techniques such as improved 



estimators ^, partial averaging higher-order ex- 

ponential splittings | [lO| , advanced reference potentials 
1 11 1, semi-classical expansions jl^], and extrapolation 
|13|. The fundamental approach is the same in all of 
these methods: the path integral is reduced to a high 
(but finite) dimensional Riemann integral, which is ap- 
proximated using either a Monte Carlo or Molecular Dy- 
namics. 

The aim of this Letter is to provide a framework for 
the formulation of a wide class of methods for the dis- 
cretization of quantum mechanical path integrals. The 
idea of approximating path integrals using a finite sub- 
set of basis functions has been suggested before in the 
literature. Davison was one of the first to consider the 
use of orthogonal function expansions in the representa- 
tion of Feynman path integrals |14|, although he did not 
explore truncating the expansion. In a related article 
on Wiener integration, Cameron proposed using a finite 
set of orthogonal basis functions, and investigated the 
convergence of Fourier (spectral) elements |15 . In this 



Letter, we do not require that the basis functions are or- 
thogonal, allowing for the direct comparison of the STA 
and FD methods. Although other authors have explored 
fundamental connections between the STA and FD |l^ ] 
methods, we are unaware of any comparison using the 
approach investigated here. In addition, our approach 
allows for the construction of new methods using gen- 
eral classes of orthogonal polynomials or finite elements. 
To illustrate the flexibility of this approach, we derive a 
new method, using compactly supported (Hermite) cu- 
bic splines (HCS), which is shown to exhibit improved 
efficiency when applied to model problems. 

To illustrate how one can use a subspace approxima- 



1 



tion to discretize the quantum density matrix in we 
start by introducing a change of variables to simphfy 
the boundary conditions and temperature dependence for 
each path integrah x (r) ^ a + {h — a) t / f3h + y [t/ fih). 
Since the admissible paths, x, satisfy the boundary con- 
ditions X (0) — a and x [(3h) = b, the reduced paths 
given by y, will satisfy Dirichlet boundary conditions, 
y (0) = y (1) — Oi independent of a, b, and (3. Intro- 
ducing this change of variables into (|l|), results in the 
following: 
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Note that the ith component of each reduced path y, de- 
noted by yi, is a real- valued function on the interval [0, 1] , 
satisfying Dirichlet boundary conditions. For the systems 
considered in this article, we also require that the deriva- 
tive of each iji is measurable (i.e., square-integrable) . 
Functions of this form are members of an infinite di- 
mensional Sobolev space defined by 5q [0, 1] = 
{w e C [0, 1] \w (0) = w (1) = and \\w\\s < oo} , where 

\\w\\l^J^w{Cf+w{0'dC 

We proceed in the following manner to discretize @): 
Consider a sequence of subspaces of increasing dimen- 
sion Vi, • • • , Vp, • • • C 5q [0, 1], where each Vp is of di- 
mension P. For convenience, let each subspace be de- 
fined as the span of a particular set of basis functions: 
Vp = spanj?/"!, • • • , ■(/'p}. Now, given a component func- 
tion Ui e Sq [0,1], we can define its projection on Vp 
uniquely by 

p 
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Using the projection yl (^) as an approximation of 
yi (^) reduces the infinite-dimensional path integral in (|^) 
to a finite-dimensional Riemann integral over the coeffi- 
cients, at k- 
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Here, we have used simplified notation for the multi- 
dimensional integral, with da = Jlfei'^o;*:,*- "^^^ con- 
stant J reflects the particular choice of variables, and 
can be readily calculated (as discussed later). The reader 
should note that (^ does not depend on the basis func- 
tions chosen to represent the approximating subspace. If 
both {^"1, ■ • • , ^p} and {^i, ■ • • , ^p} span Vp, then there 
is an invertible linear transformation (i.e., change of vari- 
ables) U such that d = Ua. 



To show in detail how subspace methods can be applied 
in practice, we consider the case of an A^-body Hamilto- 
nian system: 

1 

H = -^m,pf + V[xi,---,Xd] ■ 

Here, the coordinate and momentum operators are de- 
noted by Xi and pi respectively. For this system the 
functional $ is given by (^, which when applied to the 
projected path, x(^) (r) = a+(b-a)T//3n + y('P) (t//?^), 
results in 
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After expanding the r- integrals, introducing a change of 
variables ^ — t/ pti, and using the boundary conditions 
of each il^k , we have 



+pn yV[a+(b-a)C + y(^) (0 



where di = [ai^i ■ • • ap_i] . The "stiffness matrix", K £ 
p^PxP^ has entries given by the inner-product Kj^k = 

Substituting (^) into (||), we obtain a simplified expres- 
sion for the approximate density matrix: 



P (b, a) = exp <^ - ^ 
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For the Fourier case, one typically calculates J by requir- 
ing that the discretization be exact when applied to an 
ideal gas (i.e., V = 0) |4|j5|,[l^. Applying this same tech- 
nique to a generic subspace method, and assuming that 
K is positive definite, one can solve for J in a straight- 
forward manner: 



J = ]JVd^ 



2 7r/3r 



Before discussing particular choices for basis functions, 
we should mention that, in general, the one-dimensional 
^-integral in (^ cannot be performed analytically. This 
problem has been traditionally circumvented by using 
a discrete approximation, such as Gaussian quadrature 
1^,^. For example, one can view the primitive STA as 
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using the trapezoidal rule. If the quadrature scheme is of 
sufficiently high order its use will not reduce the asymp- 
totic rate of convergence of the overall method. An opti- 
mal scheme must be efficient, since for nonlinear iV-body 
systems evaluating V may be computationally expensive. 

As mentioned above, the real benefit of using a gen- 
eral subspace approach is the flexibility afforded through 
the choice of basis functions. By considering a general 
class of pseudo-spectral or finite-element basis functions, 
a diverse group of discretizations can be constructed. Di- 
rect comparisons can be made between basis functions of 
varying smoothness and support. However, for brevity, 
we restrict our attention in this Letter to three different 
types of basis functions: linear, spectral, and cubic ele- 
ments. Representative basis functions from each of these 
discretizations are shown in Figure |^. 



Spectral 




FIG. 1. Sample basis functions are shown above for the (a) 
linear, (b) spectral, and (c) cubic element methods. 

The traditional STA method can be constructed by 
considering polygonal paths, which can be represented by 
piecewise linear basis functions [ pT[ . For a given number 
of linear segments, we can define an approximating 

subspace Vp as the span of basis functions {V'l, ■ • • , ippji 
where each is defined by the following formula 



with 



1-H 





I)-fc), 

otherwise 



For this discretization, it is routine to show that the el- 
ements of the "stiffness" matrix, K G R^^^ , can be de- 
termined by A'^j = — l(5i_ij +2(5i.j — l(5i+i,j. In a similar 
manner, the FD method can be derived using the sub- 
space approach by considering spectral basis functions of 
the form ipk (C) 1/fcsin (fcTr^). K is diagonal for this 
basis, with entries given by Ki j = Sij. 

A new method can be constructed by approximating 
the space of paths using piecewise (Hermite) cubic splines 
(HCS) |l^]. Each spline is defined on an interval of width 
2/P, with its shape uniquely determined by its function 



value and derivative at the ends of the interval. It is 
assumed here that P is an even integer. Each piecewise 
cubic path has a continuous derivative, and is described 
by linear combinations of the basis functions 



y,hcs 
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1 < fc < P/2 
P/2<k<P 
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UG [-1,1] 
otherwise 
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One can verify that the reduced path '{^) — 
J2^k'ipk{£,) satisfies Dirichlet boundary conditions, and 
interpolates the interior grid points (2j'/P, a^) for inte- 
gers 1 < j < P/2. The derivative of the path at all 
the grid points is determined by the remaining P/2 + 1 
coefficients, ak- Due to the compact support of the ba- 
sis functions, the stiffness matrix is banded, with block 
structure: 



j^hcs 
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60 



Ki Kg 



where the blocks are given by 
Ki = 



72 -36 

-36 72 . 

. 72 -36 

-36 72 



and 
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-1 




-1 







1 4 



-3 3 . . 
0-3 . . 
. . 3 
.0-303 



Note that the blocks are not all the same size, with K3 
of dimension (P/2 — 1) x (P/2 -I- 1). The determinant of 
K^'^^ may be calculated exactly, but for most purposes it 
is enough to know that it is a constant, which will cancel 
out when (||) is used to calculate averages. 

As a numerical experiment, we apply each path inte- 
gral discretization to the problem of calculating the aver- 
age energy of a particle in a one-dimensional double- well. 
We have chosen the same double-well potential consid- 
ered in S, which is as follows: V (x) = muj'^x^/2 -t- 
A/{{x/af -1-1). The parameter values are all in atomic 
units, with u = 0.006, A = 0.009, a = 0.09, and 
TO = 1836. At low temperatures, the energy is just above 
0.006, which is below the barrier height of 0.009. 
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To measure the accuracy of each method, we compute 
the energy at a fixed temperature of T = O.lhui/k, us- 
ing MetropoUs Monte Carlo to generate the canonically 
distributed configurations. The one-dimensional line- 
integrals of the potential are approximated using Simp- 
son's rule for the FD and HCS methods, and the tradi- 
tional trapezoidal rule for the STA method. The number 
of integration nodes is set equal to the number of basis 
functions, P, resulting in the same number of potential 
evaluations for each method. For the STA method this 
results the potential is evaluated at the end points of 
each polygonal segment (consistent with its traditional 
implementation) . 

It has been previously observed that averaged quanti- 
ties (such as energy) converge at different rates, depend- 
ing on the system, reference potential, and the form of 
the estimator [pPjl^. We use a virial estimator of the 
energy j^, i? = {V{x) + xV'{x)/2), which is known to 
exhibit improved convergence properties in many prob- 
lems. The accuracy of each average is determined by 
comparing with the "exact" solution, computed by sum- 
ming over the 15 lowest energy levels as calculated with 
Numerov's method [IlOl. 
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FIG. 2. Results for the double well. 

In Fig. H the error in the computed energy is shown 
as a function of (a) the number of basis functions and 
(b) normalized CPU time. When the number of basis 
functions (or potential evaluations) is used as a measure 
of the work, we find that the FD and HCS methods are 
comparable, and both are more efficient than the STA 
method. However, when compared on the basis of CPU 
time, the HCS method is dramatically more efficient than 
both other methods. The inefficiency of the FD method 
for low-dimensional problems can be explained by consid- 
ering the work required to compute P points on the path. 
This work scales like 0{P^) for the FD method, since 
the spectral basis functions are not compactly supported. 
On the other hand, for the STA and HCS methods this 
cost scales linearly with P. Although for very high di- 
mensional problems, the cost of evaluating the potential 
should dominate, and we expect that the differences in 
computational cost would not be as pronounced. 

In summary, the problem of approximating Feynman- 
Kac path integrals can be addressed using the finite- 



dimensional subspace approach. This technique allows 
for the formulation of new methods through the choice 
of a suitable set of basis functions. In addition, tradi- 
tional methods such as the short-time approximation and 
Fourier discretization methods can be compared using 
this framework. As an example, by considering (Her- 
mite) cubic splines, a new method can be constructed 
which exhibits improved efficiency when applied to a one- 
dimensional double-well problem. 
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